﻿Module SecDifference

    Public Sub m(pic1 As PictureBox, pic2 As PictureBox, Optional KernelSize As Integer = 3, Optional ScanSize As Integer = 3)

        If KernelSize <= 0 Or ScanSize <= 0 Then
            Exit Sub
        End If

        Dim img As Matrix = New Matrix(GetPicArray(CGray(pic1.Image)))
        Dim imgH, imgW As Integer
        imgH = img.rowbound + 1 : imgW = img.colbound + 1
        Dim result1 As New Matrix(imgH, imgW)
        Dim out_h1 As Integer = imgH - KernelSize + 1
        Dim out_w1 As Integer = imgW - KernelSize + 1
        Dim temp As New Matrix(out_h1, out_w1)
        Dim center As Integer = Int(KernelSize / 2)
        For i As Integer = 0 To out_h1 - 1
            For j As Integer = 0 To out_w1 - 1
                Dim partMartix As Matrix = img.GetPartMartix(i, i + KernelSize - 1, j, j + KernelSize - 1)
                Dim v1, v2, v3, v4 As Double
                v1 = 0 : v2 = 0 : v3 = 0 : v4 = 0
                For delta As Integer = 0 To KernelSize - 2
                    v1 += (partMartix.matrix(center, delta) - partMartix.matrix(center, delta + 1)) ^ 2
                    v2 += (partMartix.matrix(delta, delta) - partMartix.matrix(delta + 1, delta + 1)) ^ 2
                    v3 += (partMartix.matrix(delta, center) - partMartix.matrix(delta + 1, center)) ^ 2
                    v4 += (partMartix.matrix(delta, KernelSize - 1 - delta) - partMartix.matrix(delta + 1, KernelSize - 2 - delta)) ^ 2
                Next
                temp.matrix(i, j) = {v1, v2, v3, v4}.Min
            Next
        Next
        temp.LetMatrix(result1.matrix, center, center)
        '将非候选点置零
        Dim num As Integer = 0
        For i As Integer = 0 To out_h1 - 1
            For j As Integer = 0 To out_w1 - 1
                If result1.matrix(i, j) <> 0 Then
                    num += 1
                End If
            Next
        Next
        Dim threshold As Double = result1.sum() / num
        For i As Integer = 0 To out_h1 - 1
            For j As Integer = 0 To out_w1 - 1
                If result1.matrix(i, j) < 10 * threshold Then
                    result1.matrix(i, j) = 0
                End If
            Next
        Next
        '获得特征点
        Dim out_h2 As Integer = imgH - ScanSize + 1
        Dim out_w2 As Integer = imgW - ScanSize + 1
        Dim result2 As New Matrix(imgH, imgW)

        For i As Integer = 0 To out_h2 - 1
            For j As Integer = 0 To out_w2 - 1
                Dim partMartix As Matrix = result1.GetPartMartix(i, i + ScanSize - 1, j, j + ScanSize - 1)
                Dim a() As Double = partMartix.Max()
                If a(0) = 0 Then
                    Continue For
                End If
                result2.matrix(i + a(1), j + a(2)) = 255
            Next
        Next
        Dim result3 As Matrix = img + result2
        result3.Normalization()
        Dim resultImg As Bitmap = ToBitmap(result3)
        pic2.Image = resultImg
    End Sub

    Public Sub f(pic1 As PictureBox, pic2 As PictureBox, Optional KernelSize As Integer = 3, Optional ScanSize As Integer = 3, Optional T_q As Double = 0.6, Optional f As Double = 0.6, Optional c As Integer = 5)
        If KernelSize <= 0 Or ScanSize <= 0 Then
            Exit Sub
        End If

        Dim img As Matrix = New Matrix(GetPicArray(CGray(pic1.Image)))
        Dim imgH, imgW As Integer
        imgH = img.rowbound + 1 : imgW = img.colbound + 1
        Dim result1 As New Matrix(imgH, imgW)
        Dim out_h1 As Integer = imgH - KernelSize + 1
        Dim out_w1 As Integer = imgW - KernelSize + 1
        Dim q As New Matrix(out_h1, out_w1)
        Dim w As New Matrix(out_h1, out_w1)
        Dim center As Integer = Int(KernelSize / 2)
        For i As Integer = 0 To out_h1 - 1
            For j As Integer = 0 To out_w1 - 1
                Dim partMartix As Matrix = img.GetPartMartix(i, i + KernelSize - 1, j, j + KernelSize - 1)

                Dim g_u2_sum As Double = 0 : Dim g_v2_sum As Double = 0 : Dim g_uv_sum As Double = 0
                For k As Integer = 0 To (KernelSize - 1) ^ 2 - 2
                    Dim row As Integer = Int(k / (KernelSize - 1))
                    Dim col As Integer = k - row * (KernelSize - 1)
                    Dim g_u, g_v As Double
                    Dim a() As Double = GetRobert(partMartix, row, col)
                    g_u = a(0) : g_v = a(1)
                    g_u2_sum += g_u ^ 2
                    g_v2_sum += g_v ^ 2
                    g_uv_sum += g_u * g_v
                Next
                Dim N As New Matrix({{g_u2_sum, g_uv_sum}, {g_uv_sum, g_v2_sum}})
                Dim trN As Double = N.Trace()
                If trN = 0 Then
                    Continue For
                End If
                Dim DetN As Double = N.Det()
                q.matrix(i, j) = 4 * DetN / (trN ^ 2)
                w.matrix(i, j) = DetN / trN
            Next
        Next

        '获取候选点
        Dim mean As Double = w.sum() / ((w.colbound + 1) * (w.rowbound + 1))
        Dim index As New Matrix(q.rowbound + 1, q.colbound + 1)
        Dim T_w As Double = f * mean
        For i As Integer = 0 To q.rowbound
            For j As Integer = 0 To q.colbound
                If q.matrix(i, j) > T_q And w.matrix(i, j) > T_w Then
                    index.matrix(i, j) = 1
                End If
            Next
        Next

        '获取特征点
        Dim out_h2 As Integer = imgH - ScanSize + 1
        Dim out_w2 As Integer = imgW - ScanSize + 1
        Dim wMatrix As New Matrix(imgH, imgW)
        Dim result As New Matrix(imgH, imgW)
        Dim t As Matrix = w * index
        t.LetMatrix(wMatrix.matrix, center, center)

        For i As Integer = 0 To out_h2 - 1
            For j As Integer = 0 To out_w2 - 1
                Dim partMartix As Matrix = wMatrix.GetPartMartix(i, i + ScanSize - 1, j, j + ScanSize - 1)
                Dim a() As Double = partMartix.Max()
                If a(0) = 0 Then
                    Continue For
                End If
                result.matrix(i + a(1), j + a(2)) = 255
            Next
        Next
        Dim result2 As Matrix = img + result
        result2.Normalization()
        Dim resultImg As Bitmap = ToBitmap(result2)
        pic2.Image = resultImg
    End Sub

    Public Sub FirDiff(pic1 As PictureBox, pic2 As PictureBox, Optional mode As String = "all")
        Dim img As Matrix = New Matrix(GetPicArray(CGray(pic1.Image)))

        Dim north_kernel As Matrix = New Matrix({{1, 1, 1}, {1, -2, 1}, {-1, -1, -1}})
        Dim east_kernel As Matrix = New Matrix({{-1, 1, 1}, {-1, -2, 1}, {-1, 1, 1}})
        Dim south_kernel As Matrix = New Matrix({{-1, -1, -1}, {1, -2, 1}, {1, 1, 1}})
        Dim west_kernel As Matrix = New Matrix({{1, 1, -1}, {1, -2, -1}, {1, 1, -1}})
        Dim north_east_kernel As Matrix = New Matrix({{1, 1, 1}, {-1, -2, 1}, {-1, -1, 1}})
        Dim south_east_kernel As Matrix = New Matrix({{-1, -1, 1}, {-1, -2, 1}, {1, 1, 1}})
        Dim south_west_kernel As Matrix = New Matrix({{1, -1, -1}, {1, -2, -1}, {1, 1, 1}})
        Dim north_west_kernel As Matrix = New Matrix({{1, 1, 1}, {1, -2, -1}, {1, -1, -1}})

        Dim kernel As Matrix
        If mode = "all" Then
            kernel = north_kernel + east_kernel + south_kernel + west_kernel + north_west_kernel + north_east_kernel + south_west_kernel + south_west_kernel
        ElseIf mode = "north"
            kernel = north_kernel
        ElseIf mode = "south"
            kernel = south_kernel
        ElseIf mode = "east"
            kernel = east_kernel
        ElseIf mode = "west"
            kernel = west_kernel
        ElseIf mode = "north_east"
            kernel = north_east_kernel
        ElseIf mode = "north_west"
            kernel = north_west_kernel
        ElseIf mode = "south_east"
            kernel = south_east_kernel
        ElseIf mode = "south_west"
            kernel = south_west_kernel
        Else
            Console.WriteLine("no such mode!!")
        End If

        Dim result As Matrix = img.Conv(kernel, 1)
        result.abs()
        result.Normalization()
        result.Binarization(12)
        Dim resultImg As Bitmap = ToBitmap(result)
        pic2.Image = resultImg
    End Sub

    Public Sub Prew(pic1 As PictureBox, pic2 As PictureBox, Optional mode As String = "all")
        Dim img As Matrix = New Matrix(GetPicArray(CGray(pic1.Image)))
        Dim vertical_kernel As Matrix = New Matrix({{1, 0, -1}, {1, 0, -1}, {1, 0, -1}})
        Dim horizontal_kernel As Matrix = New Matrix({{-1, -1, -1}, {0, 0, 0}, {1, 1, 1}})
        Dim kernel As Matrix
        If mode = "all" Then
            kernel = vertical_kernel + horizontal_kernel
        ElseIf mode = "vertical"
            kernel = vertical_kernel
        ElseIf mode = "horizontal"
            kernel = horizontal_kernel
        End If
        Dim result As Matrix = img.Conv(kernel, 1)
        result.abs()
        result.Normalization()
        result.Binarization(12)
        Dim resultImg As Bitmap = ToBitmap(result)
        pic2.Image = resultImg
    End Sub

    Public Sub Sobel(pic1 As PictureBox, pic2 As PictureBox, Optional mode As String = "all")
        Dim img As Matrix = New Matrix(GetPicArray(CGray(pic1.Image)))
        Dim vertical_kernel As Matrix = New Matrix({{1, 0, -1}, {2, 0, -2}, {1, 0, -1}})
        Dim horizontal_kernel As Matrix = New Matrix({{-1, -2, -1}, {0, 0, 0}, {1, 2, 1}})
        Dim kernel As Matrix
        If mode = "all" Then
            kernel = vertical_kernel + horizontal_kernel
        ElseIf mode = "vertical"
            kernel = vertical_kernel
        ElseIf mode = "horizontal"
            kernel = horizontal_kernel
        End If
        Dim result As Matrix = img.Conv(kernel, 1)
        result.abs()
        result.Normalization()
        result.Binarization(12)
        Dim resultImg As Bitmap = ToBitmap(result)
        pic2.Image = resultImg
    End Sub

    Public Sub SecDiff(pic1 As PictureBox, pic2 As PictureBox)
        Dim img As Matrix = New Matrix(GetPicArray(CGray(pic1.Image)))
        Dim kernel As Matrix = New Matrix({{-1, -1, -1}, {-1, 8, -1}, {-1, -1, -1}})
        Dim result As Matrix = img.Conv(kernel, 1)
        result.abs()
        result.Normalization()
        result.Binarization(12)
        Dim resultImg As Bitmap = ToBitmap(result)
        pic2.Image = resultImg
    End Sub

    Public Sub Laplace(pic1 As PictureBox, pic2 As PictureBox)
        Dim img As Matrix = New Matrix(GetPicArray(CGray(pic1.Image)))
        Dim kernel As Matrix = New Matrix({{0, -1, 0}, {-1, 4, -1}, {0, -1, 0}})
        Dim result As Matrix = img.Conv(kernel, 1)
        'Dim result As Matrix = getZeroCrossingMatrix(img.Conv(kernel, 1))
        result.abs()
        result.Normalization()
        result.Binarization(15)
        Dim resultImg As Bitmap = ToBitmap(result)
        pic2.Image = resultImg
    End Sub

    Public Sub Log(pic1 As PictureBox, pic2 As PictureBox)
        Dim newImage As Bitmap = GaussianFilter(pic1.Image)
        Dim img As Matrix = New Matrix(GetPicArray(newImage))
        Dim kernel As Matrix = New Matrix({{0, -1, 0}, {-1, 4, -1}, {0, -1, 0}})
        Dim result As Matrix = img.Conv(kernel, 1)
        result.abs()
        result.Binarization(10)
        Dim resultImg As Bitmap = ToBitmap(result)
        pic2.Image = resultImg
    End Sub

    '转为灰度图
    Public Function CGray(img As Bitmap)
        Dim h As Integer = img.Height
        Dim w As Integer = img.Width
        For i As Integer = 0 To w - 1
            For j As Integer = 0 To h - 1
                Dim c1 As Color = img.GetPixel(i, j)
                Dim cValue As Double = c1.R * 0.3 + c1.G * 0.59 + c1.B * 0.11
                Dim c2 As Color = Color.FromArgb(cValue, cValue, cValue)
                img.SetPixel(i, j, c2)
            Next
        Next
        Return img
    End Function

    '获取图像矩阵
    Public Function GetPicArray(img As Bitmap)
        Dim h As Integer = img.Height
        Dim w As Integer = img.Width
        Dim picArray(img.Width - 1, img.Height - 1) As Double
        For i As Integer = 0 To w - 1
            For j As Integer = 0 To h - 1
                Dim c1 As Color = img.GetPixel(i, j)
                picArray(i, j) = c1.R
            Next
        Next
        Return picArray
    End Function

    '矩阵转化为图片
    Public Function ToBitmap(m As Matrix)
        Dim b As New Bitmap(m.rowbound + 1, m.colbound + 1)
        Dim array(,) As Double = m.matrix
        Dim c As Color
        For i As Integer = 0 To b.Width - 1
            For j As Integer = 0 To b.Height - 1
                c = Color.FromArgb(array(i, j), array(i, j), array(i, j))
                b.SetPixel(i, j, c)
            Next
        Next
        Return b
    End Function

    '获取零交叉矩阵
    Public Function getZeroCrossingMatrix(m As Matrix)
        Dim array(,) As Double = m.matrix
        Dim temListRow As New List(Of List(Of Integer))
        Dim temListCol As New List(Of List(Of Integer))
        For i As Integer = 0 To m.rowbound - 1
            Dim temList As New List(Of Integer)
            For j As Integer = 0 To m.colbound - 1
                If array(i, j) * array(i, j + 1) <= 0 Then
                    temList.Add(j)
                    temList.Add(j + 1)
                End If
            Next
            temListRow.Add(temList)
        Next

        For i As Integer = 0 To m.colbound - 1
            Dim temList As New List(Of Integer)
            For j As Integer = 0 To m.rowbound - 1
                If array(j, i) * array(j + 1, i) <= 0 Then
                    temList.Add(j)
                    temList.Add(j + 1)
                End If
            Next
            temListCol.Add(temList)
        Next
        For i As Integer = 0 To m.rowbound - 1
            For j As Integer = 0 To m.colbound - 1
                If Not temListRow(i).Contains(j) Then
                    array(i, j) = 0
                End If
            Next
        Next
        For i As Integer = 0 To m.colbound - 1
            For j As Integer = 0 To m.rowbound - 1
                If Not temListCol(i).Contains(j) Then
                    array(j, i) = 0
                End If
            Next
        Next

        Return New Matrix(array)
    End Function

    Public Function GetRobert(m As Matrix, row As Integer, col As Integer)
        Dim array(,) As Double = m.matrix
        Dim g_u As Double = array(row + 1, col + 1) - array(row, col)
        Dim g_v As Double = array(row, col + 1) - array(row + 1, col)
        Return {g_u, g_v}
    End Function

    '高斯滤波
    Public Function GaussianFilter(b As Bitmap, Optional KernelSize As Integer = 3, Optional padding As Integer = 1, Optional sigma As Double = 1.3)
        Dim img As Matrix = New Matrix(GetPicArray(CGray(b)))
        If padding < 0 Then
            Console.WriteLine("the padding should not be a negative number!")
            Return Nothing
        End If
        Dim kernel As Matrix = New Matrix(KernelSize, KernelSize)
        For i As Integer = -padding To -padding + kernel.rowbound
            For j As Integer = -padding To -padding + kernel.colbound
                kernel.matrix(j + padding, i + padding) = Math.Exp(-(i ^ 2 + j ^ 2) / (2 * (sigma ^ 2))) / (2 * Math.PI * sigma * sigma)
            Next
        Next
        Dim kernelSum As Double = kernel.sum()
        For i As Integer = 0 To kernel.rowbound
            For j As Integer = 0 To kernel.colbound
                kernel.matrix(i, j) /= kernelSum
            Next
        Next
        Dim result As Matrix = img.Conv(kernel, padding)
        Return ToBitmap(result)
    End Function

End Module
